%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%0.01-degree precipitation producing%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Extract precipitation time series at spatial pixels
%%Note: After running the Python code"Reconstruction.py", the reconstructed IMERG-Early and SM2RAIN data for the common years and leap years are build. For each common year, a new file for date February 2 is build, by copy the data from March 1. This is designed for faster matrix calculation during the data fusion process. 
A=zeros(7011,4011,366*13,'single');
B=zeros(7011,4011,366*13,'single');
F=zeros(7011,4011,366*13,'single');

for i_year=2007:2019
    list_A=dir(['./IMERG_Early_Daily/CR_0.01/',num2str(i_year),'/','*.mat']);
    list_P=dir(['./SM2RAIN_2007_2021/CR_0.01/',num2str(i_year),'/','*.mat']);
    for i=1:366
        str= strcat ('./IMERG_Early_Daily/CR_0.01/',num2str(i_year),'/', list_A(i).name);
        data=importdata(str);
        A(:,:,i+(i_year-2007)*366)=data(start_index:end_index,:,:);

        str= strcat ('./SM2RAIN_2007_2021/CR_0.01/',num2str(i_year),'/', list_P(i).name);
        data=importdata(str);
        B(:,:,i+(i_year-2007)*366)=data(start_index:end_index,:,:);
    end
end

for i=1:7011
    A_point=A(i,:,:);
    B_point=B(i,:,:);
    save(strcat('./IMERG_Early_Daily/CR_0.01_point/',num2str(i)),'A_point');
    save(strcat('./SM2RAIN_2007_2021/CR_0.01_point/',num2str(i)),'B_point');
end


%%Double Variable Fusion
for i=1:7011
    A=importdata(strcat('./IMERG_Early_Daily/CR_0.01_point/',num2str(i),'.mat'));
    B=importdata(strcat('./SM2RAIN_2007_2021/CR_0.01_point/',num2str(i),'.mat'));
    a=Fusion(A,B);
    save(strcat('./fusion/CR_0.01_point/',num2str(i)),'x');
end


%%Producing
for year=2007:2019
    t_move=0;
    total_day=366;
    Precipitation_year=lst_365_CCI_C(:,:,366*(year-2007):366*(year-2007+1));
    if year ~= 2008 && year ~= 2012 && year ~= 2016
        total_day=365;
        Precipitation_year(:,:,60)=[];
    end
    for j=1:total_day
        day = datetime(year,01,01);
        this_day=day+caldays(t_move);
        t_move=t_move+1;
        Precipitation=Precipitation_year(:,:,j);
        this_date=datestr(this_day,'yyyymmdd');
        save(strcat('./fusion/CR_0.01_day/',num2str(year),'/',this_date),'Precipitation01');  
    end
end


fusion_folder='./fusion/CR_0.01_day/';
for year = 2007:2019
    mkdir(['./fusion/CR_0.01_day_tif/',num2str(year)])
    images = dir(strcat(fusion_folder,num2str(year),'/*.mat'));
    for dayIndex = 1:length(images)
        day = strcat(images(dayIndex).name(1:8),'_Fusion001');
        Prec=importdata(strcat(images(dayIndex).folder,'/',images(dayIndex).name));
        Prec=rot90(Prec,1);
        fileOut = strcat('./fusion/CR_0.01_day_tif/', num2str(year), '/',day, '.tif'); 
        Refference=georasterref('RasterSize',size(Prec),'Latlim',[14.95 55.05],'Lonlim',[69.95 140.05]);
        Refference.ColumnsStartFrom = 'north'; 
        geotiffwrite(fileOut,Prec,Refference); 
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%0.1-degree precipitation producing%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%Extract precipitation time series at spatial pixels
%%Note: After running the Python code"Reconstruction.py", the reconstructed IMERG-Early and SM2RAIN data for the common years and leap years are build. For each common year, a new file for date February 2 is build, by copy the data from March 1. This is designed for faster matrix calculation during the data fusion process. 
A=zeros(702,402,366*13,'single');
B=zeros(702,402,366*13,'single');
F=zeros(702,402,366*13,'single');

for i_year=2007:2019
    list_A=dir(['./IMERG_Early_Daily/CR_0.1/',num2str(i_year),'/','*.mat']);
    list_P=dir(['./SM2RAIN_2007_2021/CR_0.1/',num2str(i_year),'/','*.mat']);
    for i=1:366
        str= strcat ('./IMERG_Early_Daily/CR_0.1/',num2str(i_year),'/', list_A(i).name);
        data=importdata(str);
        A(:,:,i+(i_year-2007)*366)=data(start_index:end_index,:,:);

        str= strcat ('./SM2RAIN_2007_2021/CR_0.1/',num2str(i_year),'/', list_P(i).name);
        data=importdata(str);
        B(:,:,i+(i_year-2007)*366)=data(start_index:end_index,:,:);
    end
end

for i=1:7011
    A_point=A(i,:,:);
    B_point=B(i,:,:);
    save(strcat('./IMERG_Early_Daily/CR_0.1_point/',num2str(i)),'A_point');
    save(strcat('./SM2RAIN_2007_2021/CR_0.1_point/',num2str(i)),'B_point');
end


%%Double Variable Fusion
for i=1:7011
    A=importdata(strcat('./IMERG_Early_Daily/CR_0.1_point/',num2str(i),'.mat'));
    B=importdata(strcat('./SM2RAIN_2007_2021/CR_0.1_point/',num2str(i),'.mat'));
    a=Fusion(A,B);
    save(strcat('./fusion/CR_0.1_point/',num2str(i)),'x');
end


%%Producing
for year=2007:2019
    t_move=0;
    total_day=366;
    Precipitation_year=lst_365_CCI_C(:,:,366*(year-2007):366*(year-2007+1));
    if year ~= 2008 && year ~= 2012 && year ~= 2016
        total_day=365;
        Precipitation_year(:,:,60)=[];
    end
    for j=1:total_day
        day = datetime(year,01,01);
        this_day=day+caldays(t_move);
        t_move=t_move+1;
        Precipitation=Precipitation_year(:,:,j);
        this_date=datestr(this_day,'yyyymmdd');
        save(strcat('./fusion/CR_0.1_day/',num2str(year),'/',this_date),'Precipitation001');  
    end
end


fusion_folder='./fusion/CR_0.1_day/';
for year = 2007:2019
    mkdir(['./fusion/CR_0.1_day_tif/',num2str(year)])
    images = dir(strcat(fusion_folder,num2str(year),'/*.mat'));
    for dayIndex = 1:length(images)
        day = strcat(num2str(year),images(dayIndex).name(8:11),'_Fusion01');
        Prec=importdata(strcat(images(dayIndex).folder,'/',images(dayIndex).name));
        Prec=rot90(Prec,1);
        fileOut = strcat('./fusion/CR_0.1_day_tif/', num2str(year), '/',day, '.tif'); 
        Refference=georasterref('RasterSize',size(Prec),'Latlim',[14.95 55.05],'Lonlim',[69.95 140.05]);
        Refference.ColumnsStartFrom = 'north'; 
        geotiffwrite(fileOut,Prec,Refference); 
    end
end
